/*
 * Ray Casting CPU Orto
 *
 */

//#include "rayCastingCPU.h"
#include "Volume.h"
#include "FileManager.h"
#include "png.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

int normalizeVector3D(float * e)
{
	float magnitude = sqrtf(e[0] * e[0] + e[1] * e[1] + e[2] * e[2]);

	e[0] /= magnitude;
	e[1] /= magnitude;
	e[2] /= magnitude;

	return 0;
}

int normalizeVector4D(float * e)
{
	float magnitude = sqrtf(e[0] * e[0] + e[1] * e[1] + e[2] * e[2] + e[3] * e[3]);

	e[0] /= magnitude;
	e[1] /= magnitude;
	e[2] /= magnitude;
	e[3] /= magnitude;

	return 0;
}

float InterpolacionTri(volume_t * volume, float * p)
{
#if 0
	int x1 = (int) p[0];
	int y1 = (int) p[1];
	int z1 = (int) p[2]; 
	int x2 = x1+1;
	int y2 = y1+1;
	int z2 = z1+1;

	return 	volume->data[x1][y1][z1] * (1.0-p[0]) * (1.0-p[1]) * (1.0-p[2])	+	\
		volume->data[x2][y1][z1] * p[0]	      * (1.0-p[1]) * (1.0-p[2])	+	\
		volume->data[x1][y2][z1] * (1.0-p[0]) * p[1]	   * (1.0-p[2])	+	\
		volume->data[x1][y1][z2] * (1.0-p[0]) * (1.0-p[1]) * p[2]	+	\
		volume->data[x2][y1][z2] * p[0]       * (1.0-p[1]) * p[2]	+	\
		volume->data[x1][y2][z2] * (1.0-p[0]) * p[1]       * p[2]	+	\
		volume->data[x2][y2][z1] * p[0]       * p[1]       * (1.0-p[2]) +       \
		volume->data[x2][y2][z2] * p[0]       * p[1]       * p[2];
#endif

        int x0 = floorf(p[0]);
        int y0 = floorf(p[1]);
        int z0 = floorf(p[2]);
        int x1 = x0 + 1 < volume->dim[0] ? x0 + 1 : x0; 
        int y1 = y0 + 1 < volume->dim[1] ? y0 + 1 : y0;
        int z1 = z0 + 1 < volume->dim[2] ? z0 + 1 : z0;

//	fprintf(stdout, "(%f,%f,%f) (%d,%d,%d)  (%d,%d,%d)\n",p[0],p[1],p[2],x0,y0,z0,x1,y1,z1);

	float p000 = volume->data[x0][y0][z0];
	float p001 = volume->data[x0][y0][z1];
	float p010 = volume->data[x0][y1][z0];
	float p011 = volume->data[x0][y1][z1];
	float p100 = volume->data[x1][y0][z0];
	float p101 = volume->data[x1][y0][z1];
	float p110 = volume->data[x1][y1][z0];
	float p111 = volume->data[x1][y1][z1];

	return 	p000 * (1.0-p[0]) * (1.0-p[1]) * (1.0-p[2]) +       \
		p001 * (1.0-p[0]) * (1.0-p[1]) * p[2]       +       \
		p010 * (1.0-p[0]) * p[1]       * (1.0-p[2]) +       \
		p011 * (1.0-p[0]) * p[1]       * p[2]       +       \
		p100 * p[0]       * (1.0-p[1]) * (1.0-p[2]) +       \
		p101 * p[0]       * (1.0-p[1]) * p[2]       +       \
		p110 * p[0]       * p[1]       * (1.0-p[2]) +       \
		p111 * p[0]       * p[1]       * p[2];

}

int RayCastingOrto(volume_t * volume, unsigned char * screen, int width, int height,float max, float min)
{

	int 	i 	= 0;
	int 	j     	= 0;
	int 	z	= 0;
	int 	k 	= 0;

	float  di	= (float)volume->dim[0]/(float)height;
	float  dj	= (float)volume->dim[1]/(float)width;

	float rc = min+((max-min)/3.0);
	float gc = min+2*((max-min)/3.0);
	float bc = min+3*((max-min)/3.0);

	float * r 	= (float*) malloc(3*sizeof(float));
	float * color 	= (float*) malloc(4*sizeof(float));
	float * color0 	= (float*) malloc(4*sizeof(float));

	r[0] = 0.0;
	r[1] = 0.0;
	r[2] = 0.0;

	for (i = 0; i < height; i++)
	{
		for (j = 0; j < width; j++)
		{
			r[2] = 0;
			screen[k]    = -1;
			screen[k +1] = -1;
			screen[k +2] = -1;
			screen[k+3]  = 255;
			color[0]     = 0.0;
			color[1]     = 0.0;
			color[2]     = 0.0;
			for(z=0; z<volume->dim[2]; z++)
			{
				
				float inter = volume->data[(int)r[0]][(int)r[1]][(int)r[2]];//InterpolacionTri(volume, r);
				/*	
				color[0] += i*volume->data[(int)r[0]][(int)r[1]][(int)r[2]];
                                color[1] += j*volume->data[(int)r[0]][(int)r[1]][(int)r[2]];
                                color[2] += z*volume->data[(int)r[0]][(int)r[1]][(int)r[2]];
				*/
				if (inter > bc)
					color[2] = z*inter;
				if (inter > gc)
					color[1] = y*inter;
				if
					color[0] = x*inter;

				r[2] += 1;
			}
			screen[k]    = color[0]*255;
			screen[k +1] = color[1]*255;
			screen[k +2] = color[2]*255;
			//screen[k +3] = color[3]*255;
			fprintf(stdout, "(%u,%u,%u)   \n",screen[k],screen[k +1],screen[k +2]);
			k+=4;
			r[1] += dj;
		}
		r[0] += di;
		r[1]  = 0.0;
	}

	toPNG(screen,width,height);	

	free(r);

	return 0;
}

int main(int argc, char ** argv)
{
	int h = 150;//1280;
	int w = 150;//800;

	unsigned char * s = (unsigned char *) malloc(4*h*w*sizeof(unsigned char));;
	volume_t v;
	float pos[3] = {0.0, 0.0, 0.0};
	readHDF5(argv[1], "value", &v);
	setPositionVolume(&v, pos);

	int i = 0;
	float max = v.memory[0];
	float min = v.memory[0];
	for(i=1;i<v.dim[0]*v.dim[1]*v.dim[2];i++)
	{
		if (v.memory[i] > max)
			max = v.memory[i];
		if (v.memory[i] < min)
			min = v.memory[i];
	}
	fprintf(stdout, "Max=%f, Min=%f\n",max,min);
        float x0 = v.position[0];
        float x1 = v.position[0] + v.dim[0];
        float y0 = v.position[1];
        float y1 = v.position[1] + v.dim[1];
        float z0 = v.position[2];
        float z1 = v.position[2] + v.dim[2];
	fprintf(stdout,"Volumen (%f,%f,%f) (%f,%f,%f)\n",x0,y0,z0,x1,y1,z1);

	fprintf(stdout, "Volumen leido\n");
	
	RayCastingOrto(&v, s, w, h, max, min);

	return 0;
}
